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1.   INTRODUCTION 

Computational  fluid  dynamics  is  a  new  branch  of  physical  science 
which  has  developed  into  a  powerful  tool  with  the  advent  of  large  scale 
digital  computers.   The  complex  gas  dynamic  equations  can  now  be  solved 
numerically  to  simulate  gas  flows  under  controlled  flow  conditions. 
Computers  have  indeed  become  "numerical  wind-tunnels."  Roache  [1]  presents 
a  survey  of  previous  work  in  this  field. 

This  paper  presents  a  study  of  the  transient  flow  of  a  compressible 
viscous  gas  around  an  elliptic  cylinder  which  smoothly  accelerates  to  a 
fixed  velocity.   The  complete  time -de pendent  Navier-Stokes  equations  are 
solved  numerically  in  the  body  reference  frame  using  a  non-orthogonal  time- 
dependent  curvilinear  coordinate  system  designed  to  describe  the  body 
exactly  and  to  resolve  the  important  features  of  the  transient  flow.   The 
computation  is  performed  in  a  fixed  mesh  which  maps  into  a  grid  in  physical 
space  which  evolves  to  accommodate  the  transient  flow. 

Chapter  2  presents  the  formulation  of  the  problem  and  a  qualitative 
description  of  the  transient  flow,  and  outlines  the  development  of  a  flexible 
quasi-elliptical  time -de pendent  coordinate  system  appropriate  to  the 
problem.   Chapter  3  presents  a  description  of  the  second  order  accurate 
numerical  scheme  used  and  a  discussion  of  the  computer  implementation  of 
the  scheme.   Finally  Chapter  4  presents  numerical  results  for  Mach  number 
.169  and  local  Reynolds  number  40  for  two  cylinders,  a  nearly  circular 
cylinder  of  eccentricity  .2,  and  a  2  to  1  elliptic  cylinder  of  eccentricity 
.866.   The  results  illustrate  the  development  of  the  boundary  layer  and 


separation  and  circulation  regions  and  the  evolution  of  the  pressure  and 
flow  velocity.   The  flow  velocity  field  is  shown  in  fixed  and  accelerated 

frames  . 


2.   DESCRIPTION  OF  PROBLEM  AND  METHOD  OF  SOLUTION 

2.1.   Formulation  of  the  Problem 

Consider  an  infinite  elliptic  cylinder  with  its  major  and  minor 
axes  aligned  with  the  x  and  y  axes  centered  at  the  origin  and  surrounded 
by  a  quiescent  compressible  viscous  gas  of  infinite  volume.   The  cylinder 
starts  from  rest,  accelerates  smoothly  until  it  reaches  a  given  final 
velocity,  and  then  continues  moving  at  the  final  velocity. 

The  flow  produced  by  the  cylinder  should  have  the  following 
features.  As  the  cylinder  starts  from  rest,  a  weak  wave  is  produced  which 
propagates  at  the  speed  of  sound.  As  the  ellipse  accelerates,  mass 
accumulates  in  front  of  the  body  and  cavitation  occurs  in  the  rear.   When 
the  final  velocity  is  reached,  the  mass  accumulated  in  front  of  the  body 
produces  a  compression  wave.   Similarly,  the  void  in  the  rear  produces 
an  expansion  wave.   This  phenomenon  is  very  similar  to  the  one  dimensional 
piston  problem[2].  As  these  waves  propagate,  they  interact  in  a  complex 
way,  and  eventually  a  steady  state  is  reached  near  the  body.   Our 
objective  is  to  simulate  this  flow  by  solving  the  complete  time -de pendent 
Navier-Stokes  equations  in  the  body  reference  frame. 

The  smooth  acceleration  of  the  cylinder  permits  the  simulation 
of  the  transient  gas  flow  as  well  as  the  steady  flow  around  the  body. 
In  previous  studies  of  exterior  flows  described  by  the  Navier-Stokes 
equations,  the  body  was  either  started  impulsively,  or  the  steady  inviscid 
flow  was  used  as  the  initial  flow.   These  studies  were  aimed  at  only  the 
steady  flow  solution  because  the  transient  flow  produced  by  these  initial 


conditions  is  physically  meaningless.  An  exception  is  the  study  of  the 
transient  inviscid  flow  around  a  circular  cylinder  by  Moretti[3]. 

2.2.   Adaptive  Coordinate  System 

Cartesian  coordinates  are  best  suited  for  rectangular  bodies 
because  the  grid  can  be  designed  so  that  the  body  coincides  with  grid 
lines.   For  curved  bodies,  however,  a  rational  choice  of  body  nodes  leads 
generally  to  an  irregular  Cartesian  grid,  and  the  body  is  only  approximated 
by  a  piecewise  linear  curve  joining  the  body  nodes.   The  accuracy  of 
this  approximation  can  be  increased  by  increasing  the  number  of  body  nodes, 
but  the  price  paid  is  a  large  number  of  unnecessary  nodes  and  further 
irregularity  in  the  grid.   The  situation  is  further  complicated  by  the 
existence  of  a  boundary  layer  which  requires  additional  nodes  near  the 
body  for  resolution. 

These  problems  can  be  resolved  through  the  use  of  a  curvilinear 
coordinate  system  in  which  the  body  coincides  with  grid  lines.   The 
natural  system  for  this  study  is  the  orthogonal  elliptical  cylindrical 
coordinate  system  (u,v)  in  the  accelerated  body  reference  frame.   This 
coordinate  system  has  other  advantages  besides  exact  body  representation. 
A  grid  in  (u,v)  with  uniform  u  and  v  spacing  maps  into  a  grid  in 
physical  space  whose  nodes  cluster  near  the  body.   The  concentration  of 
points  is  highest  where  the  curvature  is  greatest.   These  two  characteristics 
of  the  coordinate  system  combine  to  place  more  nodes  in  physical  space 
where  the  gradients  in  the  flow  variables  are  greatest.  A  typical  grid 
corresponding  to  a  uniform  (u,v)  grid  for  an  ellipse  of  eccentricity  .866 


is  shown  in  Figure  1. 

This  elliptical  coordinate  system  would  be  ideal  for  monitoring 
small  perturbations  from  the  steady  state  solution.   However,  the  use  of 
this  coordinate  system  for  the  transient  flow  studied  here  leads  to 
unresolvable  problems.   First,  as  waves  move  away  from  the  body,  they 
move  into  the  area  where  the  nodes  are  widely  spaced,  and  hence  they 
cannot  be  resolved  adequately.   Second,  at  the  fixed  outer  boundary, 
signals  arriving  from  the  body  cannot  be  handled  properly.   The  compression 
and  expansion  waves  exhibit  steep  gradients  which  must  be  resolved  by  a 
locally  fine  grid,  and  the  signals  arriving  at  the  outer  boundary  must  be 
handled  carefully  to  prevent  the  creation  of  spurious  waves  which  can 
destroy  the  calculation.  Adding  nodes  in  the  neighborhood  of  the  waves 
and  at  the  outer  boundary  would  be  one  solution,  but  it  is  impractical 
because  of  computer  space  and  time  limitations. 

A  solution  is  to  introduce  a  time-dependent  non-orthogonal  quasi- 
elliptical  coordinate  system  designed  so  that  the  nodes  of  a  fixed 
rectangular  grid  in  this  system  map  into  points  in  physical  space  which 
move  to  accommodate  the  transient  flow.   First  u  is  mapped  into  a  new 
variable  Q   by  a  bilinear  transformation  so  that  the  interval  along  any 
appropriate  v  hyperbola  from  the  body  to  the  outer  boundary  is  mapped  into 
[0,1].   Then  Q   is  mapped  into  the  final  variable  u*  using  a  hyperbolic 
transformation.   The  outer  boundary  initially  is  an  ellipse  concentric 
with  the  body  and  of  the  same  eccentricity  but  slightly  larger.   In  the 
fixed  frame,  the  outer  boundary  expands  uniformly  at  the  speed  of  sound, 
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asymptotically  approaching  a  circle  centered  at  the  origin.   Hence  the 
signals  created  by  the  moving  body  remain  within  the  computational  space. 
As  the  calculation  proceeds  the  transformation  parameters  adjust  the 
mapping  so  that  the  fixed  grid  nodes  map  into  points  in  physical  space 
clustered  about  the  steep  gradients  in  the  flow. 

2.3.   Coordinate  System 

The  transformation  from  the  fixed  frame  (x,y)  to  the  accelerating 
body  frame  (x*,y*)  is  given  by 


x*  =  x 


r  w  dt  , 

d   x 


t 

y*  =  y  -  J*  w   dt   , 
o 

where  w  and  w  are  the  x  and  y  velocities  of  the  body  in  the  fixed  frame, 
x      y 

By  varying  w  and  w  ,  the  body  can  move  in  any  direction  at  any  angle 
attack  and  at  any  speed.   The  conservative  form  of  the  Navier-Sokes  equations 
in  the  accelerating  frame  with  density,  momentum,  and  energy  as  the 
primary  variables  are  presented  in  Appendix  A. 

The  transformation  from  the  coordinates  (x*,y*)  to  the  elliptical 
coordinates  (u,v)  is  given  by 

x*  =  c  cosh(u)  cos(v)   , 
y*  =  c  sinh(u)  sin(v) 

The  u  and  v  level  curves  are  confocal  ellipses  and  hyperbolas  with  foci 


(-c,0)  and  (c,0).  An  ellipse  of  eccentricity  e  is  represented  by  the  u 

level  curve  for  u  ,  where 
o 

uq  =  ln[(l+(l-eV)/e]   . 

For   this   ellipse   the    focal    length  is   c,    the   semi -major  axis   is   c   cosh(u   ), 

and  the  semi -minor  axis  is  c  sinh(u  ).   The  coordinate  u  ranges  from 

u   to  °°  and  v  ranges  from  0  to  2tt.   Hence  this  transformation  maps  the 
o 

exterior  of  the  ellipse  into  a  semi-infinite  strip.   The  region  of 
interest,  however,  is  a  rectangle  which  lengthens  in  the  u  direction  with 
time.   The  Navier-Stokes  equations  of  Appendix  A  transformed  to  elliptical 
coordinates  are  presented  in  Appendix  B.   In  the  transformed  equations, 
the  terms  due  to  the  accelerated  reference  frame  are  retained  in  cartesian 
form  for  convenience . 

The  transformation  from  u  to  Q   is  given  by  the  bilinear  mapping 

C  =  (bu+c)/(u+a) 


where 


a  =  [(2uo-u1)uoo-u1uo]/(2u1-uo-uJ   , 


b  =  (Ul+a)/[2(Ul-uo)]   , 


c  =  -  u  b 
o 


This  transformation  maps  the  inner  boundary  u  ,  an  intermediate  point  u, , 
and  the  outer  boundary  u  into  0,  1/2,  and  1.   Hence  it  maps  the  region 


of  interest   into  a    fixed  rectangle. 

The   final   transformation  is   given  by 


*  -  I  +  i_  i    ri+(2C-l)tanh(g(t)/2)-| 
u      "  2        lot       Ll-(2C-l)tanh(a(t)/2)J      ' 


v*  =  v      , 


t*  =  t 


This  transformation  preserves  the  fixed  rectangular  computational  space, 
but  is  compresses  the  region  in  Q   near  1  and  0  and  stretches  the  region 
near  1/2.   The  computational  space  is  covered  by  a  fixed  uniform 
rectangular  grid  which  maps  into  a  grid  in  physical  space  whose  nodes  are 
clustered  in  the  boundary  layer,  particularly  in  the  regions  of  greatest 
body  curvature,  and  near  the  outer  boundary.   The  parameters  u.  and  a 
are  used  to  control  the  distribution  of  the  nodes  in  physical  space. 

The  Navier-Stokes  equations  are  solved  on  this  uniform  grid. 
All  factors  in  the  equations  other  than  the  derivatives  of  the  flow 
variables  can  be  easily  evaluated  analytically.   The  partial  derivatives 
of  the  flow  variables,  however,  must  be  approximated  by  finite  differences 
The  first  partial  derivatives  of  a  function  f  are  written  in  terms  of 
u*  and  v*  by  applying  the  chain  rules 


f  -  f  ,.  u*r  C    , 

u     u*    fa   u 


f   =  f   +  f  ,  u*  £ 
v    v*    u*   C  v 


ft  =  ft*  +  fu*(U*C  Ct  +  u*a  at' 
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The  second  partial  derivatives  are  transformed  by  applying  these  equations 
a  second  time.   Then  the  derivatives  of  flow  variables  with  respect  to 
u*,  v*,  and  t*  are  approximated  by  finite  differences  on  the  uniform  grid. 

2.4.   Node  Distribution 

In  the  body  reference  frame,  the  compression  and  expansion  waves 
move  at  different  speeds.   Hence  the  outer  boundary  at  u,,,  is  a  function 
of  v*  and  t*.   In  the  fixed  Cartesian  frame,  the  outer  boundary  is  represented 

by 

\xu+c  t  /  +  \y,+c  t  /    1      ' 
bo       b   o 

where  c   is  the  speed  of  sound  in  the  gas  at  rest,  and  x,  and  y,  are  the 
o  b      b 

semi -major  and  semi -minor  axes  of  the  boundary.   Transforming  this  representa- 
tion to  the  accelerated  frame  moving  only  in  the  x  direction  and  scaling 
distances  with  respect  to  c  cosh(u  )  yields  the  following  equation  for  u^, 

u  (v*,t*)  =  ln[A+(A2-l)2]   , 


where 


CO 


2     i 
A  =  (-C+(C  -4BD)2)/(2B)   , 

2         2  2 

B  -  (cos  (v)E+  sin  (v)F)/(cosh  (u  )E)   , 

o 

C  =  (2  x,  cos(  v))/cosh(u  )   , 

G  O 

D  .=  x,2  -  (F  sin2  (v))/(E  cosh2  (u  ))  -  F   , 
a  o 

e  =  (yb+cot)2  , 

F  =  (x,  +  c  t)2 

D     O 
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Here  x,  is  the  nondimensiona 1  distance  the  ellipse  has  moved. 

The  parameters  <*(t*)  and  u  (v*,t*)  control  the  node  distribution 
The  parameter  u..  is  chosen  so  that 

cosh(u  )  -  cosh(u  ) 

.___ 1_  =  \ 

cosh(u^)  -  cosh(u  ) 

where  X  is  a  constant  between  0  and  1.   As  X   increases  nodes  are  pulled 
away  from  the  outer  boundary;  as  \   decreases,  nodes  are  pulled  away 
from  the  body. 

Solving  this  equation  for  u,(v*,t*)  yields 


Ul(v*,t*)  =  ln[A+(A2-  1)2]  , 


where 


A  =  (l-\)cosh(um)  +  A,  cosh(u  ) 

o 

Since  the  outer  boundary,  u  ,  expands  while  the  total  number 
of  nodes  remains  fixed,  the  nodes  tend  to  move  away  from  the  body.   However, 
the  grid  in  physical  space  must  be  fine  enough  to  resolve  the  boundary 
layer.   Hence  a  maximum  length  6  for  the  distance  in  physical  space 
between  the  body  and  a  first  u*  node  from  the  body  at  some  value  of  v* 
is  specified  for  each  simulation.   This  length  6,  obviously  is  inversely 
proportional  to  the  Reynolds  number. 

The  parameter  a(t*)  is  used  to  control  the  spacing  of  the  first 
row  of  u*  nodes  from  the  body.   The  initial  value  of  a   is  chosen  so  that 
the  maximum  distance,  6,  is  attained  at  a  given  time  t  .   For  times 


12 


smaller  than  t   or  is  held  constant  and  the  spacing  increases  until  it 
equals  6  at  time  t. .   Thereafter,  a   is  chosen  to  keep  the  spacing  equal 
to  6. 

By  carefully  choosing  X  and  6,  the  node  distribution  can  be 
tailor  to  the  problem  needs.   For  example,  in  a  low  Reynolds  number  flow 
containing  strong  waves,  X  should  be  small  to  concentrate  nodes  near 
the  waves.   In  a  high  Reynolds  number  flow  containing  weak  waves,  X 
should  be  large  to  concentrate  nodes  near  the  body  and  6  should  be  small 
enough  to  permit  adequate  resolution  of  the  boundary  layer. 

A  typical  node  distribution  in  the  (u*,v*,t*)  space  for  some 
time,  t  >  0,  is  shown  in  Figure  2. 
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3.      NUMERICAL  METHOD  AND  IMPLEMENTATION 

The  Navier-Stokes  equations  in  the  (u*,v*,t*)  coordinate  system 
can  be  written  in  the  form 


f  .  -  G(f  .,f  „  , ,f  ,,f  ,  , ,f  ,  ..) 
t*      u*  u*u*  v*  v*v*  u-"v* 


where 


T 
f  =  {p,m  ,m  ,E)    , 

p  =  density   , 

m  =  u  component  of  momentum   , 
u 

m  =  v  component  of  momentum   , 
v 

E  =  stagnation  internal  energy 

These  equations  were  solved  numerically  using  MacCormack's  predictor- 
corrector  method [4].   This  scheme  was  chosen  because  it  is  second  order 
accurate,  is  easy  to  implement,  and  has  worked  well  on  viscous  flow 
problems  . 

Because  of  computer  space  and  time  limitations,  the  angle  of 
attack  is  set  equal  to  0  to  exploit  the  symmetry  of  the  flow  about  the  x-axis 
The  body  velocity  w  is  described  by  the  cubic  spline  with  knots  0,  tf/2, 
and  t   satisfying 

w  (0)  =  0   , 
x 

wx(t£/2)  =  v£/2   , 
wx(t)  -  v£,   t  >  t£-  . 
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This  spline  is  continuous,  has  two  continuous  derivatives,  and  is  anti- 
symmetric about  (v_/2,t_/2). 

The  computational  space  in  (u*,v*)  is  the  rectangle  [0,1]  x  [0,tt] 

The  grid  used  is  {u*,v*}  where 

i  j 


u*  -  iAu*,   1  ■  0,  1,  ...,  N  , 
i 

v*  =  jAv*,   j  -  0,  1 N   , 

j 

where  Au*  =  l/M  and  Av*  =  tt/N.   The  time  is  partitioned  by  the  grid  {t*} 

n 


where 


e*-0 

t*   =  t*  +  At*,   n  =  0,  1,  . . .   . 
n+1    n     n 

The  value  of  any  function  f(u*,v*,t*)  at  the  point  (u*,v*)  at  time  t*  is 

i   j  n 

denoted  by  f .  . . 

3.1.   Exterior  Flow  Field 

The  exterior  flow  field  is  defined  by  {u*,v*},  i  =  l,  •••  ,  M,  and 

i      j 

j=0,    ...    ,   N.      MacCormack's    predictor-corrector  method   first   predicts   the 

values   of   the    flow  variables   at   time   t    , ,    using 

n+1 

_n+l  n  n 

f.     .   =  f.     ,   +  At  G.     .,      i=l,    ...    ,   M-l  and  j=0,    ...    ,   N      , 
i,J  i.J  i,J 

where 

n  n  n2"2n  n 

G.     ,   s  G(A     f.     .,    Av*f.     . ,    6   ^f.     .,    6      fn      ,   A  ^  ^f        )       , 
i,j  u*  i,j'  i,j'      u*  i,j'      v*  i,j'      u*v*  i,2 
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A    .f .     .    =    (f,  .,     ,    -    f.     .)/Au*      , 
u*  i,j  i+l,  J  t,J 

A  *f4     •    =    (f-     ,xi    "    f-     -)/Av*      . 
v*  i,j  i,j+l         1,J 

62,f.     .   *   (f4    .        -   2f.     .   +  f.  .,    .)/Au*2      , 
u*  i,j  l-l, j  i,j  i+l, j 

62   f.     .   =   (f.     ,    ,    -  2f,     .   +  f.     .,.)/Av*2      , 

A    ,.  ,f.     •   =   (f-  L1       «_,    -    f-^-.        -   f-       ...    +  f-     .)/Au*Av*      . 
U*V*    l,j  1+1,  j+1  1+1,  j  1,J+1  i,j' 

All   the    terms   in  G  are   evaluated  at   time   t*  and   the    predicted  values   of   the 

n 

flow  variables  at  the  outer  nodes  are  set  to  the  freestream  conditions  at 


time  t*  . .. 
n+1 


Then  the  final  values  of  the  flow  variables  at  t  ,.  are  obtained 

n+1 


using  the  following  correcter  equation, 

n+1     n     _n+l     _n+l  , 

f.  .  =  (f.  .  +  f.  .  +  AtG.  ,)/2,  i=l,  ...  ,  M-l  and  1=0,  ...  ,  N, 
i.J     i-.J    i,J      i>J 

where 

-n+1        -n+i     _n+l   2  -n+1   .2  ~n+l  „    ^+lx 
G.  ,  =  G(V  f.  .,  V  f.  .,  6  .f.  .,  6  J..     .,  V    f.  .)   , 
i,j     v  u*  i,j'   v*  i,j'   u*  i,j'   v«  i,j'   u*v*  i,j 


V  f   ,  =  (f.    -  f.  .  .)/Au*   , 
u*  l,j      i,j     i-l, J 


V  f.  .  =  (f.  .  -  f.  . 

v*  1,J       l,j     1,J 


f.    J/Av*  > 


62  f.  .  =  (f.  .,  .  -  2f.  .  +  f  -   .)/Au*2   , 
u*  i,j     i+l, j     i,j    i-l, j 
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62    f      .    =   (f.     ...    -  2f,     .   +  f.     .    .)/Av*2      , 
v*  i,j  i,j+l  i,j  i,J-l 


?  +_*,£.     ,   s   (f.        ~   f-     ,    i    ~   f.    n    •   +  f.    ,        J/Au*Av*      . 
u*v*  i,j  i,j  i,j"l         i-l, J  1-1, J-l' 

—  "k 

All   the   terms   in  G  are   calculated  at   time   t    .,   and  the   final  values   of  the 

n+1 

flow  variables  at  the  outer  nodes  are  set  to  freestream  conditions  at  time 


'n+1 


Both   the   predictor  and  the  corrector   require   the   flow  variables 


at  v*  =   -   Av*  and  v      .    =  rr  +  Av*.      These  values   are   obtained  by  symmetry 

from  the  values   at  v*  and  v„   ,    as    follows 

i  N-l 

Pi,j    =  Pl,j+2 


(m  ).     .    =    (m   ),     .,„ 


(m   ) .     .    =   -    (m   ) .        „ 
vi,j  v  i,j+2 


E.     .    =  E.     ,  ,0 
i,J  i.j+2 

for   j    =   -    1  and  N-l. 

3.2.      Flow  at   the   Body 

A  separate  calculation  is   required  to  compute   the   flow  variables 
at   the   body.      The   two  major  conditions   imposed  on  the   flow  at   the  body  are 
(1)   no-slip  at   the   surface   in   the  accelerated   frame,   and    (2)    the  wall  is 
adiabatic . 
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The  no-slip  or  zero  particle  velocity  condition  is  implemented 

by  setting  m  and  m  to  zero  at  the  wall.   For  an  adiabatic  wall,  the 
u      v 

internal  energy  of  the  gas,  e,  at  the  surface  has  zero  gradient  normal  to 
the  surface. 

u  body 
Now  the  energy  of  a  perfect  gas  is 

E  =  pe  +  p  V-V/2   . 
Since  the  velocity  vector,  V,  at  the  wall  is  zero,  it  follows  that 

(E  ),  ,   =  (p  e).  ,   =  (p  E/p).  . 
u  body     u  body     u    body 

To  be  consistent  with  MacCormack's  second  order  scheme,  all 
derivatives  must  be  approximated  by  a  second  order  accurate  difference 
formula  at  the  wall.   Hence  the  following  equation  is  used 

(f  )   .  =  (l/2Au*)(-3f   .  +4f.  .  -  f,  ,)(u*C  )   .   • 
u  o,j  o,j     l,j    2, j    q     u  o,j 

Substituting  this  expression  for  E  and  p   in  the  equation  already  given 
and  solving  for  E   .  yields 

E   .  =  p   .(4E.  .  -  E0  .)/(4p1  .  -  p0  .)   . 
o,J    o,jv   l,j    2,j     l,j   r2,j 

In   the    predictor   p       ,    can  be   calculated  using   the   continuity 

°,J 

equation  and  the  predicted  energy,  E   . ,  is  obtained  from  the  equation  for 

°  >  J 

E   .  above  using  predicted  values  at  time  t* , ,  on  the  right  hand  side. 
o,j  °  n+1 
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The  boundary  calculations  in  the  corrector  present  a  problem 
since  nodes  (u*  ,v*)  are  needed  but  do  not  exist.   The  required  derivatives 
in  the  continuity  equation  must  be  approximated  by  a  second  order  accurate 
difference  formula  consistent  with  MacCormack's  corrector  equation.   The 
following  approximation  can  be  used. 

(f  )   .  =  (l/Au*)(-2f   .  +  3f .  .  -  f ,  ,)(u*C  )   . 
u  o,j  o,j     l,j    2,j   q    u  o,j 

Using  this  result  in  the  continuity  equation  yields  the  values 

of  p   ..   The  energy,  E   . ,  is  then  obtained  from  the  equation  for  E 

o,J  o,j  o,j 

above  using  corrected  values  at  time  t*  .  on  the  right  hand  side. 

A  second  method  recommended  by  Roache[5]  for  computing  the  density 
at  the  boundary  was  also  used.   This  method,  which  employs  a  staggered 
mesh,  provided  results  very  close  to  the  results  of  the  first  method.   The 
method  required  a  large  amount  of  geometric  computation  because  of  the 
time  dependent  grid  and  moreover  required  the  solution  of  a  tridiagonal 
system  of  linear  equations  for  the  density  values.   Since  this  method  for 
all  its  added  complexity  did  not  produce  significantly  different  results, 
it  was  discarded  in  favor  of  the  first  scheme. 

3.3.   Computer  Simulation 

The  complexity  of  the  coordinate  system,  Navier-Stokes  equations, 
and  numerical  scheme,  the  small  time  steps  required  for  stability,  and 
the  relatively  large  number  of  nodes  required  for  adequate  resolution, 
all  combine  to  require  a  very  large  amount  of  computer  time  for  each 


20 


simulation.   Another  important  factor  is  the  amount  of  internal  storage 
available.   Since  many  time  dependent  variables  in  the  predictor  are 
identical  to  those  of  the  previous  corrector,  these  variables  can  be 
computed  at  each  node  once  and  used  both  in  the  corrector  at  t  and  the 
predictor  at  t    if  storage  is  available.   Unfortunately,  storage  was  at 
a  premium  on  the  machines  used,  and  many  quantities  had  to  be  recomputed. 
This  added  to  the  already  large  computing  times  for  each  simulation. 

Initially  the  outer  boundary  is  a  finite  distance  from  the  body. 
Hence  calculations  are  performed  only  at  the  nodes  within  the  envelope 
of  the  wave  created  by  the  initial  motion  of  the  cylinder.   At  the  remaining 
nodes,  the  flow  variables  are  assigned  their  freestream  values.   This 
procedure  saves  a  large  amount  of  computer  time  during  early  times,  and 
works  better  than  the  alternative  procedure  of  adding  nodes  using 
interpolation. 

The  numerical  scheme  was  implemented  in  a  Fortran  program.   Two 
versions  of  the  program  were  written:   one  for  the  PDP-10  timesharing 
system,  and  one  for  the  IBM  360/75  batch  system.   The  basic  programs  were 
identical,  but  the  PDP-10  program  contains  interactive  facilities  which 
allow  the  user  to  review  results  dynamically,  set  checkpoints,  and  adjust 
the  time  step  for  stability.   Both  versions  automatically  adjust  the  time 
step  to  maintain  stability.   If  instability  occurs  the  time  step  is 
decreased  and  the  program  restarts  using  the  data  at  the  previous  stable 
checkpoint.   Checkpoints  for  program  restarts  are  created  by  writing  the 
current  values  of  all  time  dependent  variables  into  a  disk  file.   This  pro- 
cedure minimizes  the  effect  of  a  machine  failure  or  numerical  instability. 
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The  parameters  required  for  a  simulation  are  the  reference 
Reynolds  number,  eccentricity,  grid  parameters  M  and  N,  node  distribution 
parameters  X,  6,  x,  ,  and  t  ,  acceleration  parameters  t.  and  Vf,  and  the 
initial  time  step  At*.  When  a  run  is  restarted  using  a  checkpoint  of  a 
previous  run,  At*  may  be  reset.   The  results  are  printed  on  a  line 
printer  and  plotted,  if  desired,  on  a  CALCOMP  plotter  or  a  Tektronix 
4010-1  display. 
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4.   NUMERICAL  RESULTS 

Calculations  were  performed  for  two  cylinders,  a  nearly  circular 
cylinder  of  eccentricity  .2  and  a  2  to  1  cylinder  of  eccentricity  .866. 
In  both  calculations  Au*  is  1/31,  Av*  is  tt/37,  \   is  .68,  6  at  v*  =  tt  is 
.008,  the  Prandel  number  is  .75,  y   is  1.4,  the  reference  Reynolds 
number  is  100,  the  final  Mach  number  is  .169  attained  at  time  .5,  and 
the  local  Reynolds  number, 

Re.  =  2c  cosh(u  )p  w  (t)/u   , 

I  O   O   X 

is  40  at  Mach  number  of  .169. 

In  both  calculations  the  initial  time  step  is  .0001.  As  the 
minimum  spacing  between  nodes  increases  to  6,  the  time  step  increases  to 
.001  for  the  first  cylinder  and  .00095  for  the  second.   Because  of  the 
small  time  steps  used,  each  calculation  required  a  large  number  of  steps 
and  consequently  a  large  amount  of  time.   The  first  required  15.5  hours 
on  thelBM  360/75  to  reach  time  11.4  in  13,362  steps,  while  the  second 
required  18  hours  to  reach  time  13.46  in  15,380  steps'. 

As  the  body  accelerates,  mass  accumulates  in  front  of  the  body 
and  cavitation  occurs  behind  the  body.   This  phenomenon  is  evident  in 
the  3-dimensional  pressure  plots  in  (u*,v*)  shown  in  Figures  4  and  5, 
where  the  pressure  is  high  in  the  front  and  low  behind.   These  two 
figures  also  illustrate  the  effects  of  viscosity  since  they  are  from  runs 
differing  only  in  the  reference  Reynolds  number.   The  viscosity  in  Figure  5 
is  10  times  that  in  Figure  4.   Consequently  the  smearing  of  the  compression 
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Figure  3.  Description  of  3-dimensian  pressure  field  plots  in 
(u*,v*)  space. 
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Figure  4.   3 -Dimensional  pressure  plot,  eccentricity  is  .866, 
Re  is  1000,  Mach  number  is  .034,  and  at  time  .14. 


Figure  5.   3-Dimensional  pressure  plot,  eccentricity  is  .866, 
Re  is  100,  Mach  number  is  .104,  and  at  time  .29. 
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Figure  6.   3-Dimensiona 1  pressure  plot,  eccentricity  is  .866, 
Re  is  100,  Mach  number  is  .169,  and  at  time  2.95. 


Figure  7.   3-Dimensiona 1  pressure  plot,  eccentricity  is  .866, 
Re  is  100,  Mach  number  is  .169,  and  at  time  10.25. 
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Figure  8.   3 -Dimensional  pressure  plot,  eccentricity  is  .866, 
Re  is  100,  Mach  number  is  .169,  and  at  time  13.46. 


Figure  9.   3-Dimensiona 1  pressure  plot,  eccentricity  is  .2, 

Re  is  100,  Mach  number  is  .169,  and  at  time  10.60, 
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and  expansion  waves  is  greater  in  Figure  5. 

When  the  final  velocity  is  attained,  compression  and  expansion 
waves  move  out  from  the  front  and  rear  of  the  body  at  the  speed  of  sound, 
Vv.   As  these  waves  interact,  they  create  complex  waves  within  the  gas. 
These  waves  can  be  seen  in  Figures  6  and  7.   Eventually  a  steady  state 
is  reached  near  the  body.   Figures  8  and  9  show  the  steady  state  pressure 
distribution  near  the  body  as  well  as  the  compression  and  expansion 
waves  far  from  the  body  for  both  cylinders. 

The  pressure  coefficient  is  defined  as 

Cp  =  2(p-l)/„x2   , 

where  p  is  the  dimensionless  pressure  and  w  is  the  dimensionless  velocity 

of  the  body.   The  evolution  of  C  at  the  surface  of  the  bodies  is  shown 

P 

in  Figures  10  through  13.   Initially  the  pressure  is  high  in  front  and 
low  behind,  but  as  the  compression  and  expansion  waves  travel  around  the 
body  and  interact,  the  pressure  decreases  in  front  and  increases  behind 
and  the  pressure  coefficient  reaches  the  steady  state  shown  in  Figure  11 
at  time  10.60  for  the  near  circle  and  in  Figure  13  at  time  13.46  for  the 
2  to  1  ellipse . 

The  density  variation  throughout  the  flow  field  is  7%  for  the 
near  circle  and  5%  for  the  2  to  1  ellipse.   The  temperature  variation 
at  the  body  is  270  for  the  near  circle  and  1%  for  the  2  to  1  ellipse.   The 
variations  for  the  near  circle  are  larger  than  those  for  the  2  to  1 
ellipse  because  the  blunter  near  circle  produces  stronger  compression 
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Figure  10.   Pressure  coefficient  around  body  of  eccentricity  .2. 
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Figure  11,   Pressure  coefficients  around  body  of  eccentricity  .2 
for  various  times.   Time  10.60  is  steady  state  curve. 
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Figure  12.   Pressure  coefficient  around  body  of  eccentricity 
.866. 
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and  expansion  waves.   This  difference  in  relative  strengths  of  the  waves  is 
evident  in  the  3-dimensiona 1  pressure  plots  in  Figures  8  and  9. 

Figures  14  through  22  show  the  velocity  field  development  in 
the  body  reference  frame.   They  illustrate  how  the  boundary  layer 
evolves  and  how  separation  and  circulation  eventually  develop.   Figures  18 
and  22  show  that  separation  occurs  further  downstream  for  a  slender  body 
than  for  a  blunt  body.   Let  s  be  the  arc  length  along  the  body  from  the 

front  stagnation  point  and  let  I   be  half  the  body  perimeter.   Figures  18 

s 
and  22  show  separation  occurring  at  y  =  .63  for  the  near  circle  and  at 

s 

j   =  .84  for  the  2  to  1  ellipse. 

Figures  23  through  29  show  the  velocity  field  development  in  the 
fixed  frame,   They  illustrate  how  the  accelerating  body  affects  the  gas 
particles.   The  particles  at  the  surface  move  in  the  same  direction  and  at 
the  same  speed  as  the  body.   The  particles  in  front  of  the  body  are  pushed 
out  and  around  the  cylinder,  while  those  behind  are  pulled  in  toward  the 
body  to  fill  the  void.  A  circulation  region  develops  near  v  =  tt/2  during 
acceleration,  and  then  moves  slowly  downstream.   Figures  23  to  29  illustrate 
the  development  of  these  regions. 

These  results  represent  nearly  incompressible  isothermal  flow. 
In  principle  this  method  can  also  be  applied  to  high  subsonic  Mach  number 
flows  where  compressible  and  thermal  effects  are  more  significant.   However, 
finer  meshes  would  be  required  to  resolve  the  steeper  gradients  of  such 
flows,  and  smaller  time  steps  would  be  required  for  stability.   This  might 
raise  the  computer  storage  and  processor  time  to  prohibitively  high  levels 
on  many  machines . 
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Figure  14.  Initial  boundary  layer  growth  at  time  .29 
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Figure  15.   Velocity  field  development  at  time  3.04 
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Figure   16.      Velocity  field  development  at   time   10.25 
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Figure    17.      Circulation  and  separation  at   time    13.46 
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Figure  18 


•    0f  circulation  region  at  time  13.46, 
Enlarged  view  of  circus 
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Figure  19.  Initial  boundary  layer  growth  at  time  .38, 


Figure  20.  Velocity  field  development  at  time  2.69 
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Figure  21.   Circulation  and  separation  at  time  10.60, 
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Figure  22.  Enlarged  view  of  circulation  region  at  time 
10.60. 
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Figure  23.      Initial  gas    particle  movement  at   time    .29 
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Figure  24.   Gas  particle  movement  at  time  3.04. 


Figure  25.   Gas  particle  movement  at  time  10.25 
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Figure  26.   Circulation  region  of  gas  particles  has  moved 
downstream  at  time  13.46. 
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Figure  27.     Initial  gas    particle  movement  at   time    .38 
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Figure  28.      Gas   particle  movement  at  time  2.69. 
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Figure  29.   Circulation  region  of  gas  particles  has  moved  down- 
stream at  time  10.60. 
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APPENDIX  A 
Nayier-Stokes   Equations,   Accelerated  Reference   Frame 

The    two   dimensional  Navier-Stokes   equations    for  an  ideal  com- 
pressible viscous   gas   in  an  accelerated  Cartesian  reference   frame    (x*,y*) 
are   presented  here. 

The    transformation   from  the   fixed   frame    (x,y)    to   the  accelerated 
frame    (x*,y*)   is   given  by 

t 

x*   =  x    -    f  w     dt      , 
•'      x 
o 

y*  =  y  -  J  w    dt     , 

o      y 

where  w  and  w  are  the  x  and  y  velocities  of  the  accelerated  frame  relative 
x      y 

to  the  fixed  frame . 

The  primary  dependent  variables  chosen  are  the  density  p,  the  x* 

and  y*  components  of  the  momentum  m  ^  and  m  .,  and  the  total  energy  E. 

The  momentum  and  total  energy  are  defined  as  follows: 


V  =  pyy*    ' 

E  =  P(e+(V^+V^)/2)   , 

where  V  .  and  V  .  are  the  x*  and  y*  components  of  the  gas  velocity,  and  e 
is  the  internal  energy.   Other  dependent  variables  are  the  velocity  V, 
temperature  T,  stress  t,  pressure  P  and  heat  flux  of  q. 
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All  dimensional  quantities  will  be  nondimensionalized  with  respect 
to  the  following  characteristic  quantities. 

L  =  length   , 

o  ■  velocity   , 


where 


Ia/y  = 


time   , 
c 
o 


P   =  density 


PC2 

o  o  , 
=  pressure  and  energy   , 


2 
c 

-—  =  temperature   , 


c  =  specific  heat  at  constant  volume   , 
v 


c  =  specific  heat  at  constant  pressure   , 
P 

cp 
Y  =  ratio  of  specific  heats,  —  , 

v 

R  -   gas  constant,  c   -  c    , 

P    v 


c  =  speed  of  sound  in  reference  gas 
o 


The  scaling  yields  the  following  nondimensional  quantities 
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p  c  L 
Re  ■  reference  Reynolds  number 


Pr  =  Prandtl  number  =  -r*   , 

M  ■  Mach  number  =  —   , 

c 
o 

where  k  is  the  thermal  conductivity  and  \i   is  the  viscosity. 

The  Navier-Stokes  equations  in  dimensionless  form  with  the  stars 
dropped  for  convenience  are: 
Continuity 

-    dm    oni 

&  +  — *  +  — £  =  0 

ot   dx    By       ' 

x -momentum 


dm    -.  ^  -, 

TT  +  =T  (m  V  )  +  f-  (m  V  )  =  5-  [t  -P  ] 
ot    ox   x  x    oy   x  y    ox   xx 


St      dw 
xy       j 

+  o7^"  p  dT 


y -momentum 


Smv   *  a  * 

TT-  +   I"  (m  V  )  +  f-  (m  V  )  =  ~  [T   -P] 

ot    ox  s  y  x    ay   y  y    oy   yy 


oTxy     dwy 

+  Sx~"  ■  P  dT 
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Total  energy  §£  +  ^  (EV,)   +  ^   (EVy)   -  -  ^  -  Jl 


+  §-    [V     (T        "  ?)   +   V     T        ] 

ox       x     xx  y    xy 


+  j—  [V      T         +  V    (T       -p)] 

oy       x     yx         y     yy 


dw  dw 

x  _y 

x  dt  y  dt 


where 


T§§      Re  Lo?        Vd5    +  oT7/3J     » 

T§^]  "  Re   LoTl     +  d§  J      ' 


Y  ST 

q?  "         (Y-DPrRe    &5      ' 


for   5  =  x,y  and    M  =  y,x. 


Note   that  w     and  w     are   relative   to  the   fixed   frame. 
x  y 
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APPENDIX  B 
Navier-Stokes  Equation  in  Elliptic  Coordinates 

This  appendix  presents  the  Navier-Stokes  equations  of  Appendix  A 
transformed  to  elliptic  cylindrical  coordinates.   The  transformation  from 
the  Cartesian  coordinates  (x,y)  to  the  elliptic  coordinates  (u,v)  is  given  by 

x  =  c  cosh(u)  cos(v)   , 
y  =  c  sinh(u)  sin(v) 

The  notation  and  variable  definitions  are  identical  to  those  of  Appendix  A. 
However,  the  reference  length  L  of  Appendix  A  is  now  assumed  to  be  the 
semi -major  axis  of  the  cylinder,  c  cosh(u  ),  and  the  momentum  is  now  resolved 
into  u  and  v  components.   For  convenience,  the  x  and  y  components  of  the 
body  velocity  w  and  w  are  retained  in  Cartesian  form. 


The  transformed  Navier-Stokes  equations  [6]  are 


Continuity 


ft +  ^ [§;*%> +  t;<h°v]  =  0  • 


u -momentum 


Bm 


■ST1  =  ^r  1?-   [h(T       -m     V     -T        +  m     V)] 
ot         ,  2    Lgu  uu        u     u         vv         v     v 


+  2  f-   [h(T       -m     V  )]    +  h 
ov  uv         u     V 


r-   [T        -    P  -  m     V   ] 
.ou    W  V     V 


a  11  dwx  dwy 

§-  [T        -m     V  ]  \Y  +  Apl-  +  BP  jr-      , 
ov       uv         u     v  JJ  dt  dt 


v -momentum 
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Total  energy 


3m 

a*1  -  -T  if"   ft(T        -m     V     -t       +m     V  )] 
ot         ,2    ^-ov    w         v     v         uu  u     u 


5T  ~  .2  lav  [hl 


+  2  t-   [h(T        -  m 
ou  vu         v 


V   )]    +  h||-   [T        -  P   -  m     v  1 
u  Ldv       uu  u     uj 


ou        vu  V 


m     V   ]])  - 
V     u  JJ 


dw  dw 


If  ■  h  ih  <h<VTuU 


-  P 


E)    +  V      T 

V     vu 


Vi 


+  f~   [h(V   (t        -   P   - 
OV  V  w 


E)    +  V     T 

u     uv 


V1 


where 


(      dwx  dwy\       (     dwy  dwx\ 

+  AVm     TT~  +  m     ~a7~)  +  BVm     7Z m     77/     ' 

\  u  dt  v  dt    /  \  u  dt  v  dt    ' 


2,    .     .      .2,    ,,2 


(sinh    (u)    +  sin    (v))2/cosh(u   )      , 


h 

A  =  cos(v)    8inh(u)/(h   cosh(u   ))      , 

B  =  sin(v)    cosh(u)/(h  cosh(u   ))      , 


SV.       V-, 


.     .  2_  { i  Zl  +  3  Ml 

§1       Re   lh   a?         u2   ST| 


77  rlf  <hV  +  Iff  <hV]}     ■ 
3n 
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-  1     \ 


m-  fc  t&  (Yh)  +  tn  (Vh)) 


ST 


q§   "   "    (Y-l)PrRe  h  d§      ' 


for   5  ■  u,v  and  T]  ■  v,u. 
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